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Abstract. We study the phase transition between survival and extinction in an 
epidemic process with long-range interactions and immunization. This model can 
be viewed as the well-known general epidemic process (GEP) in which nearest- 
neighbor interactions are replaced by Levy flights over distances r which are distributed 
as P{r) ~ r~'^~°'. By extensive numerical simulations wc confirm previous field- 
theoretical results obtained by Janssen et al. [Eur. Phys. J. B 7, 137 (1999)]. 



1. Introduction 

In statistical physics models for epidemic spreading continue to attract attention because 
they exhibit non-equilibrium phase transitions with universal properties. The best 
known example is the universality class of directed percolation (DP) [1-4], which 
is observed in lattice models for epidemic spreading with nearest-neighbor infection 
and spontaneous recovery. Very recently, this type of phase transition was observed 
experimentally for the first time by Takeuchi et al. [5]. 

In the past two decades various generalized models of epidemic spreading have been 
proposed. One possible generalization is to introduce a local memory in order to mimic 
the effect of immunization [6]. This leads to the so-called general epidemic process 
(GEP), which has been studied both numerically and by field-theoretic methods [7-9]. 
Another possibility is to consider long-range infections which are usually modelled by 
so-called Levy flights (for a recent review see [10]). In such models the infection is 
transported over randomly chosen distances r, which are power- law distributed as 

P(r)~r-'^"", (1) 



where d is the spatial dimension and a > is an exponent controlling the characteristic 
shape of the distribution. By studying DP with Levy flights it turns out that critical 
exponents vary continuously in a certain range of a. 
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Let us briefly recall some of the main results for the GEP. Starting point 
is the Langevin equation for directed percolation, i.e., epidemic spreading without 
immunization. 

dt(l){x, t) = a0(f , t) - 60^ (f, t) + DV^(j){x, t) + ^{x, t). (2) 

Here is the density of infected individuals, a is a parameter controlling the spreading 
rate, — reflects to lowest order the constraint that individuals cannot be doubly 
infected, and D'V'^(f) is a diffusion term accounting for nearest-neighbor interactions. The 
noise ^ describes density fluctuations caused by the stochastic nature of the spreading 
process and is deflned by the density- dependent correlations 

{(f){x,t)(f){x',t')) = r(l){x,t)S'^{x-x')S{t-t'). (3) 

As flrst suggested by Cardy [7], the effect of immunization can be implemented by 
adding a term of the form 

X(f){x,t) exp (^-^ dt' (j){x,t')^ . (4) 

In this term the integral sums up the past activity at position x between the initial 
condition t — and the actual time t. If this integrated activity is still small 
the exponential function is essentially equal to 1, marking a non-immune individual. 
However, when the integrated activity exceeds a certain threshold of the order 1/w, the 
exponential function 'switches' to zero, representing an immunized individual. Omitting 
arguments x, t the modified Langevin equation reads 

dtcj) = a(p- h(p^ + DVV + C + A exp ^ dt' 0(f , t')^ . (5) 

As can be seen, the exponential function is coupled to the field 0(a;, t) so that it effectively 
modifies the spreading rate a. This means that non-immune sites are infected with rate 
a -|- A, while immune sites are infected with rate a. This allows one to control the rates 
for the first and all subsequent infections of an individual separately. The additional 
term modifies the critical behavior of the transition, leading to a universality class which 
is different from DP. This so-called GEP-class comprises all models which are defined 
in the spirit of this Langevin equation. 

Using the Janssen-de Dominies formalism by introducing a response field 0(a;, t) 
and integrating out the Gaussian noise, one is led to a field-theoretic action [7, 8] 

j d'^x j dt 4>(^-dt + a + iPV')*/) + - (6) 

+ ^4>(f) exp (^j^ dt' (f){x, t' 
which is expected to be valid in arbitrary dimensions. 
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Table 1. Critical exponents of the standard general epidemic process. Empty entries 
are not available. 



Apart from the interpretation as an epidemic process with immunization, the GEP 
also plays an important role as a dynamical process that produces isotropic (undirected) 
percolation clusters in d dimensions. More specifically, whenever the process terminates, 
it leaves behind a certain cluster of immune sites. This cluster can be shown to be an 
ordinary percolation cluster [11] for which the critical exponents are already known (in 
two dimensions even exactly) and can be expressed in terms of two standard exponents 
/3 and u = i'±. In the GEP, time as an additional degree of freedom induces another 
exponent, the so-called dynamical exponent z. Knowing (3, u, and z, the exponents for 
the survival probability P(t) ~ and the average number of particles N(t) ~ can 
be expressed as 

S = —, e = -- — -l for ordinary GEP. (7) 

Z/|, Z 

wehre i/y = zux. The upper critical dimension of the GEP is cic = 6 and the exponents 
are summarized in Table [H 

In this paper we consider the GEP generalized in such a way that the nearest- 
neighbor interactions are replaced by non-local ones according to a Levy distribution ([T]). 
For large values of cr, where the Levy fiights become short-ranged, one expects to recover 
the standard GEP. On the other hand, for sufficiently small values of a the interactions 
are so long-ranged that any spatial structure is wiped out, hence one expects mean- field 
behavior. In between there is an intermediate regime described by a non-trivial field 
theory [14], where critical exponents vary continuously with a. 

The paper is organized as follows. In the next section we recall how the problem 
can be formulated as a field-theoretic action that involves a fractional derivative. Then 
in Sect. [3] we analyze the scaling behavior of this action, determine the range where 
perturbative methods can be applied, discuss two scaling relations which are conjectured 
to hold to all orders of perturbation theory, and quote the one-loop results derived by 
Janssen et al. [14]. Finally we verify these findings by numerical simulations in Sect. HI 
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Figure 1. Typical snapshot of clusters of immune sites for different values of a. 



2. Field-theoretic action 



A Levy flight in d dimensions is an isotropic random displacement x — > x + s according 
to a probability distribution which for large distances decays algebraically as 



P{s) 



~ \s 



-d—a 



(8) 



Here the exponent cr > is a control parameter which determines the asymptotic power- 
law characteristics of long-range flights. In order to normalize this distribution, a lower 
cutoff of this power law for small s is needed, for example in form of a minimal flight 
distance. 



Acting on a density fleld t), a Levy flight generates the time evolution 

Introducing the non-local linear operator V*^ which acts on a function f{x) as 

1 



(9) 



V-/(f) 



d'^slsl"'^-" 



with the (7-dependent normalization constant 

7r'^/2r(-|) 



f{x + s)-f{x) 



(10) 



this equation of motion can be written to leading order as 



(11) 



(12) 



which has the same structure as the ordinary diffusion equation apart from the fact 
that the Laplacian is replaced by V'^. In the literature the operator V'^ is known as a 
fractional derivative because it has certain algebraic properties that generalize those of 
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ordinary derivatives. For example, for < o" < 2 the action of the operator V'^ on a 
plane wave amounts to bringing down a prefactor of the form 



V" e 



ik-x 



We 



ik-x 



(13) 



As a field-theoretic approach to GEP with Levy flights, one may simply replace the 
Laplacian in the action ([6]) by a fractional derivative. However, past experience in 
the study of directed percolation with Levy flights [14-17] shows that the fractional 
derivative should be added as a new term without discarding the Laplacian. The reason 
is that even if the Laplacian was not included in the bare action, it would be generated 
under renormalization group. Therefore, the field-theoretic action describing the long- 
ranged GEP is given by 



d'^x / dt 



b~ 



(14) 



^ 2 



exp 



dt' 0(f, t') 



3. Scaling analysis 



3.1. Scaling dimensions 



The major problem technical problem when investigating this action lies in the structure 
of the exponential function. In order to apply ordinary field-theoretic methods, one 
needs to expand the exponential function as a power series. However, such an expansion 
is only meaningful if the relevance of the terms under renormalization group decrease 
with increasing order. 

As usual one introduces the dimensions 

[x] = [k]-\ [t] = [k]-\ [4>] = [k]\ [^] = [kr (15) 

with respect to a momentum scale k, where z = is the dynamical exponent and 

X, X ci-re the field exponents. As usual, the fields and (p stand for the creation and 
removal of active sites. Since each removed particle leaves an immune site behind, 
the scaling properties of the generated cluster of immune sites are described by the 
exponent x- Therefore, one can identify x with the exponent (3 of isotropic percolation 
through 



(16) 
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3.2. Expansion of the exponential function 



Obviously, the integrand in the last term of the action has the dimension 



dt' 0(f , t') 



[k] 



X-2 



(17) 



This means that the exponential function can be expanded as a power series only if 
X > z since in this case the terms of the series are decreasingly relevant. Under this 
condition the zeroth term of the series, ^(fxp, can be absorbed by a shift a a + X/2 
so that the non-trivial term to be retained in the action is the first-order term of the 
exponential |00 jjdt'0(x, t') which has the dimension [fc]^"*"^^"^. For ^ > this term 
is more relevant than the cubic term — in Eq. so that the latter one can be 
discarded. Hence for x > z > the action takes the form 



Srr 



d^x / dt 



-dt + a + DV' + D\/''](j) + - 



+ 



Xw ~ 



dt'Mx, t') 



For a = 2 this expression reduces to the action of the conventional GEP. 

By dimensional analysis of the field-theoretic action (fT8|) one obtains the mean field 
(tree level) exponents 



1 



1/a 



1 



X 



2(T, 



-MP 



a. 



(19) 



which are expected to be valid above the upper critical dimension 
dr = 3a . 



(20) 



3.3. Exact scaling relations 

The exponents X ^i-iid z are related by two exact scaling relations [14]. These relations 
can be explained as follows. 

The first scaling relation is caused by a rapidity reversal symmetry of the action. 
As can be verified easily, the action of the standard GEP is invariant under the 
replacement [7,9] 

|0(f,t) ^ 0(f,-t). (21) 

Obviously this rapidity reversal symmetry holds also in the case of long-range 
interactions a ^ 2. Moreover, it is an exact symmetry which holds to all orders of 
perturbation theory. Comparing the scaling dimensions on both sides this immediately 
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implies the exact scaling relatioi||| 

X = X-z. (22) 

As long as the field exponent x is positive (which is indeed the case here), the condition 
X> z is fulfilled, justifying the expansion of the exponential term in Eq. flTSl) . 

The second scaling relation can be made plausible as follows. The term with Levy 
operator D(j)V^(f) is bilinear in the fields. Therefore, it modifies the free theory, whereas 
it does not change the interactions described by the cubic terms. In other words, the 
Levy term modifies the free propagator of the theory, but it has no infiuence on the 
structure of the Feynman diagrams. This means that all loop integrals have exactly the 
same form as in the standard GEP, the only difference being that the free propagator in 
momentum space = a + iu + Dk'^ has to be replaced by = a + iu + Dk'^ + Dlkl'^ . 
In so far the situation is completely analogous to the case of directed percolation with 
long-range interactions. 

In the case of Levy-DP it was observed that the corrections caused by the loop 
integrals can be expressed as a power series in the momentum k. As in any Taylor 
expansion, this series involves only integral powers, the leading order being k'^. This 
means that the long-ranged Levy operator produces loop corrections that can be 
absorbed in the coefficient of the (short-range) Laplacian. In other words, the Levy term 
is not renormalized, but instead renormalizes the short-range diffusion coefficient D. 
This implies that the Levy term with scaling dimension ^k]^~^^~^'^~'^~^ must be invariant 
under rescaling, hence one obtains the scaling relation 

X + x + cr-d-z = (23) 

This so-called hyperscaling relation holds to all orders of perturbation theory below the 
upper critical dimension d < dc- Combining it with Eq. (|22|) one obtains 

2x = rf-a, (24) 

or equivalently 

/3 = ^(rf-a). (25) 
The hyperscaling relation implies that the exponents 

0='-^±^-l (26) 

Z/|| Z/|| 

are related by 

7^ = ^- (27) 
9 + 1 2a ^ ^ 



X Note that this scaUng relation differs from the one for directed percolation, where the well-known 
rapidity reversal symmetry (f>{x, t) <-> (f>{x, —t) forces x &n.d x to be identical. 
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3.4- Crossover from long-ranged to ordinary GEP 

The scaling relations are particularly interesting at the crossover from the regime, where 
the Levy operator is relevant, to the ordinary GEP, where the short-ranged Laplacian 
plays the dominant role. This crossover takes place at a particular value of a denoted 
as as- Assuming that the crossover is smooth in the sense that the exponents do not 
change discontinuously, as can be computed by plugging the known GEP exponents 
into the hyperscaling relation (l23l) and solving it for a: 

'^s = d-^. (28) 

For example, in = 2 dimensions the threshold for a, above which the short-range GEP 
is recovered, is given by 

43 

a. = — ^ 1.792 (rf = 2) . (29) 

At a first glance this result seems to be counterintuitive as it contradicts the naive 
expectation that the fractional derivative V'^ should reduce to the ordinary Laplacian 
at the threshold a = 2. However, studies of directed percolation and the Ising model 
with Levy-type interactions have shown that this naive argument does not necessarily 
hold in an interacting theory. Moreover, it is surprising that in 2D the crossover takes 
place at a threshold cr^ < 2, while in previous studies of directed percolation as was 
always found to be larger than 2. 



3.5. One-loop results 

According to Janssen et al. [14], the critical exponents computed in d = dc — e = "ia — e 
dimensions to one-loop order are given by 



/? =l-^ + 0(e^). 


(30) 


^--^4>^(^^)' 


(31) 


,,^i+,;+o(e^ 


(32) 




(33) 




(34) 


' -16>^^^^)- 


(35) 



Usually a field-theoretic e-expansion can be verified only approximately by numerical 
simulations: since the dimension d is an integer number, it is in most cases impossible to 
study small values of e -C 1. However, in the present model the upper critical dimension 
dc = 3(7 is controlled by the parameter a and can be chosen in such a way that e 
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Figure 2. Critical percolation threshold Ac of the long-range general epidemic process 
as a function of a (see text). 



becomes small. This allows one to test the one-loop results directly through Monte 
Carlo simulations, as we show in the next section. 



4. Numerical analysis 



4-1- Strategy 

For a simulation in d dimensions, let us introduce the parameter /i by 

^ ■= a-d/3. (36) 

For /i = the (i-dimensional system is at the upper critical dimension, while for a small 
/i > it is slightly below dc with e = d^ — d = 3fi. Substituting this parameter the 
exponents 6 and 6 are given to one-loop order by 

Hence, plotting these exponents as functions of a, one expects the following behavior: 

• In the mean field regime a < d/3 the exponents are constant and given by their 
mean field values. 

• In the interval d/3 < a < cTs the exponents vary continuously. The slope of the 
tangent at the left edge is just one-loop correction in Eq. (1371) . 

• Finally, for a > the Levy operator becomes irrelevant and the exponents of the 
standard GEP are recovered. 



4-2. Details of the simulation 

The model studied in our simulations is defined on a two-dimensional square lattice 
of sites that can be in the states 'healthy' (H), 'contaminated' (C) and 'immune' (/). 
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The model evolves by the following microscopic transitions: On the one hand, each 
contaminated site becomes immune at a rate A 

On the other hand, a contaminated site may infect another site, which is randomly 
selected by a Levy flight at distance r: 

C ^-"y- H C C 

These transitions were implemented by random-sequential updates on the two- 
dimensional grid with periodic boundary conditions and lateral length L. An integer 
variable s{x,y) is attached to each lattice site {x,y) and stored in a two-dimensional 
array. s{x, y) = means that the site is in state H, s{x, y) = 1 stands for a site in state 
C and s{x, y) = 2 denotes a site in state I. In addition, we kept all infected sites in 
a dynamically generated list to accelerate the numerical calculations. The dynamics of 
the simulation is generated by a repetition of the following steps: 

1) Randomly choose one out of the list of infected sites. 

2) With the probability A the site is immunized by the transition C —>■ I, 
followed by a jump to step 4. 

3) Otherwise, a second site is selected by a Levy flight according to the a- 
dependent distribution ([1]) (see below). If this target cell is in the 'healthy' 
state, it will be infected by the transition from state H to C. 

4) If the number of infected sites is zero, the simulation stops, otherwise 
increase the system time t by the reciprocal total number of infected sites 
and go back to step 1. 

Hence the model is controlled by two parameters, namely, the Levy exponent a > and 
the critical parameter A G [0, 1]. 

Upon infection, the isotropic Levy flight used to locate the second site is distributed 
by a power-law P(r) ~ j--"^-^ where r = a/ Aa;^ + Ay"^ is the absolute value of the 
distance between two sites. In two dimensions, the radial distance r and the equally 
distributed angle can be generated by setting r := and := 1'nz2^ where 

-2^1,2 ^ [0, 1] are random numbers. Note that the radial distribution has a lower cutoff 
so that the minimal flight distance is 1. The actual coordinates are then computed by 
rounding the cartesian coordinates (Ax, Ay) = (r cos0, r sin 0) to integers. 

In order to avoid the time-consuming use of trigonometric functions, one may 
accelerate simulations significantly by doing the following: first one repeatedly generates 
a random vector in a square with real- valued coordinates Axs,A|/s G [0,1] until 
= Ax^ + Ay^ < 1. This procedure generates an isotropically distributed random 
vector on the unit disk and with a squared radius equally distributed between and 1. 
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Then this radial distribution has to be corrected in such a way that the desired Levy 
distribution ([1]) is obtained. This can be done by setting 

Ax:=^-r = Ax, - (Ax,' + Ay,')-'^'-'^" 

Ay:=^.r = Ay,. (Ax/ + Ay,')-'^'-'^'' (38) 

Levy flights exceeding the system size are treated by applying periodic boundary 
conditions. As discussed in Ref. [15], alternative methods where such flights are 
dismissed or truncated are plagued by much stronger finite-size effects. 



4-3. Numerical results 



In order to estimate the critical exponents of anomalous DP we used dynamical 
simulations at criticality, where a critical cluster is grown from a single active seed 
(just as in Fig. [1]). In our simulations, the initial condition was a L x L - lattice with 
healthy sites, in which we placed four neighboring contaminated sites forming a seed. 
Averaging over many independent realizations we measured the survival probability P{t) 
and the actual number of contaminated sites N{t). At criticality, these quantities are 
expected to scale as 

P{t) oc , N{t) oc , 

where 6 = P/u^ and 9 = {dh'± — 2[3)/i'\\ — 1. Since deviations from criticality lead to a 
curvature of P(t) and N{t) in double logarithmic plots, we held a G [0.2, 7.0] constant 
for different values of A, searching for a straight line in order to get a precise estimate 
of the percolation threshold A for different values of a (see Fig. [2]). For small a, finite 
size effects even occur at small simulation times, which can be minimized by increasing 
the lateral system size L. For a = 0.2 and L = 14000, as an example, finite size effects 
occurred for T > 2000. For bigger a, we simulated up to T = 10000 and averaged over 
5000 runs. 

As can be seen from Fig. [3] one observes three different regimes as expected (see 
discussion after equation (|371)): 

• cr < 2/3: mean-field behavior (interactions range is effectively infinite). 

• 2/3 < a < as'- intermediate phase with continuously varying exponents whose 
values, to one-loop approximation, are given by (l37|l . 

• a > as'- ordinary short-range GEP with exponents given in Tab. [H where the Levy 
Operator becomes irrelevant. 

Our numerical results are shown in Fig. [3], where the expected values of the mean-field 
exponents and the short-range GEP exponents are marked by horizontal lines. The 
field-theoretic one-loop predictions according to (137|) are shown as tilted straight lines, 
which in principle should be tagent to the numerical curve near the crossover. The 
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coincidence is far from being perfect, probably because of strong finite-size effects when 
cr is small. Nevertheless the estimates seem to follow the predicted behavior and, as far 
as we can jugde, confirm the field-theoretic results. 
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5. Conclusions 



In this paper we have studied the general epidemic process (GEP) with long-range 
interactions and immunization. In our model infections are transported by Levy flights 
over large distances according to a power-law distribution P(r) ~ f~'^~'^ depending on 
a parameter a > 0. After being infected once, individuals become perfectly immune. 

The main result of our work is that we fully conflrm, through numerical analysis, 
the fleld-theoretic predictions of Janssen et ai, namely that one should observe three 
different regimes depending on the value of a. For a < 2/3, the system is so strongly 
mixed that mean-field behavior is expected at the transition. On the other hand, 
for a above a characteristic threshold cTg ^ 1.792 the usual short-range Laplacian 

dominates over the Levy operator which, in field-theoretic parlance, becomes 
irrelevant. In this regime the system becomes effectively short-ranged and one recovers 
GEP exponents. In between, one has an intermediate regime with continuously varying 
scaling exponents. To one- loop order these exponents are given by ( |3711 . in fair agreement 
with the numerical simulations for small values of o" > 2/3. 

The (T-dependent simulation is interesting in that it allows us to move continu- 
ously from the mean-field regime over the long-range phase to the short-range regime 
without changing the space dimension d. At the crossover from the mean field to the 
LGEP phase, it is possible to verify the field theoretic one-loop prediction quantita- 
tively. The second crossover from the LGEP phase to the short-range GEP regime is 
also interesting in so far as it takes place, for = 2, at cr^ ~ 1.792 and not at a = 2, 
as one would naively expect as a threshold where replaces the ordinary Laplacian V^. 
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